1. Simulation tests

To determine how long to run our simulations, we ran a series of tests to evaluate when population size, fitness, and genetic structure reached equilibrium. We ran each unique set of simulation parameters (i.e., varying population size, selection strength, migration, spatial autocorrelation, and correlation between environmental layers) for 10,000 timesteps. We simulated 1,000 neutral loci in the test simulations (in contrast to 10,000 in the full simulations) to decrease computational time; the smaller number of loci should not notably influence the results of any of our equilibrium tests. We extracted population size every 1 timestep, mean fitness every 10 timesteps, and genetic data every 1,000 timesteps (the genetic data files are large and computationally expensive to process). Based on the results of our tests, described in detail below, the latest our simulations reached equilibrium was ~4,000-5,000 timesteps. Therefore, we decided to run our simulations for 6,000 timesteps to ensure they had enough time at equilibrium.

2. Testing for equilibrium in population size

We tested for equilibrium in population size using an Augmented Dickey-Fuller test for stationarity. We tested for stationarity in windows of 1,000 time points (e.g., 0-1,000, 1,000-2,000, etc.) for each unique set of simulation parameters. We calculated the percentage of simulations that were at stationarity within those windows based on an alpha level of 0.05 (alternate hypothesis: stationarity). Population size acheived stationarity by 1,000 timesteps.

3. Testing for equilibrium in fitness

3.1 Stationarity test

We tested for equilibrium in mean fitness using an Augmented Dickey-Fuller test for stationarity. As before, we tested for stationarity in windows of 1,000 time points and calculated the percentage of simulations that were at stationarity within those windows based on an alpha level of 0.05. The majority of the simulations were at stationarity by 2,000 timesteps, however the percent of simulations at stationarity didn’t reach a maximum until the 9,000-10,000 timestep window (A). This is likely due to very slight variation in the variance of fitness that is unlikely to affect our inferences about sampling strategy. We confirmed that all of the simulations had reached stationarity in at least one of the stationarity windows by 4,000 timesteps (B).

3.2 Plateau test

To determine the timepoint at which fitness plateaued, we fit a linear plateau model to the mean fitness data. We used this model to deterimine the junction point at which fitness flattened out (i.e., the timepoint at which the slope of the model changes to zero). The latest junction point was at ~2,000 timesteps.

4. Testing for genetic equilibrium

We tested for equilibrium in genetic structure using two tests: (1) a test for equilibrium in genetic structure using Mantel tests and (2) a test for equilibrium of genotype-environment associations based on genotype-environment correlations.

4.1 Genetic structure test

To determine when genetic structure reached equilibrium, we calculated the correlation between (1) genetic distance and geographic distance and (2) genetic distance and environmental distance using Mantel’s r. We calculated these correlations every 1,000 timesteps for each unique set of simulation parameters. There are not enough points to run an Augmented Dickey-Fuller test for stationarity, so we visually inspected the data to determine when the correlations reached equilibrium. Each line on the plot represents a unique simulation; we colored the lines based on the migration and spatial autocorrelation parameters of each simulation because they are expected to have strong effects on neutral and adaptive genetic structure, however these lines also include variation in population size, selection strength, and environmental correlation. Overall, the Mantel correlations are stable by about 4,000-5,000 timesteps.

4.2 Genotype-environment association test

To determine when genotype-environment associations reached equilibrium, we calculated genotype-environment correlations for both adaptive and neutral loci. We calculated correlations every 1,000 timesteps for each unique set of simulation parameters. We calculated correlations two seperate ways: (1) including fixed loci and counting them as a correlation of 0 and (2) excluding fixed loci. In the plots including fixed loci, correlation for neutral loci decreases overtime as more loci become fixed due to the loss of genetic diversity to drift overtime in the absence of mutation (B). This decrease is unlikely to affect our analyses because we simulated 10,000 loci, so many will still be variable, even at later timepoints. The adaptive loci plots (A) show greater variation over time because they are averages across just 8 loci (compared to 1,000 neutral loci for the neutral loci plots). The plots are otherwise the same as described above; each line on the plot represents a unique simulation, colored by migration and spatial autocorrelation parameters. Overall, the correlations for the adaptive and neutral loci appear to be relatively stable by about 4,000 timesteps.